home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / bivnor.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  3KB  |  140 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. #define twopi 6.283195307179587
  15. #define con (twopi / 2.0) * 10.0e-10
  16.  
  17. double bivnor(ah, ak, r)
  18.      double ah, ak, r;
  19. {
  20.   /*
  21.     based on alg 4628 comm. acm oct 73
  22.     gives the probability that a bivariate normal exceeds (ah,ak).
  23.     gh and gk are .5 times the right tail areas of ah, ak under a n(0,1)
  24.  
  25.     Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand.
  26.   */
  27.   double a2, ap, b, cn, conex, ex, g2, gh, gk, gw, h2, h4, rr, s1, s2, 
  28.     sgn, sn, sp, sqr, t, temp, w2, wh, wk;
  29.   int is;
  30.  
  31.   temp = -ah;
  32.   normbase(&temp, &gh);
  33.   gh = gh / 2.0;
  34.   temp = -ak;
  35.   normbase(&temp, &gk);
  36.   gk = gk / 2.0;
  37.  
  38.   b = 0;
  39.  
  40.   if (r==0)
  41.     b = 4*gh*gk;
  42.   else {
  43.     rr = 1-r*r;
  44.     if (rr<0)
  45.       return(-1.0);
  46.     if (rr!=0) {
  47.       sqr = sqrt(rr);
  48.       if (ah!=0) {
  49.     b = gh;
  50.     if (ah*ak<0)
  51.       b = b-.5;
  52.     else if (ah*ak==0)
  53.       goto label10;
  54.       }
  55.       else if (ak==0) {
  56.     b = atan(r/sqr)/twopi+.25;
  57.     goto label50;
  58.       }
  59.       b = b+gk;
  60.       if (ah==0)
  61.     goto label20;
  62.     label10:  
  63.       wh = -ah;
  64.       wk = (ak/ah-r)/sqr;
  65.       gw = 2*gh;
  66.       is = -1;
  67.       goto label30;
  68.     label20:  
  69.       do {
  70.     wh = -ak;
  71.     wk = (ah/ak-r)/sqr;
  72.     gw = 2*gk;
  73.     is = 1;
  74.       label30:  
  75.     sgn = -1;
  76. /*    t = 0;  not used JKL */
  77.     if (wk!=0) {
  78.       if (fabs(wk)>=1)
  79.         if (fabs(wk)==1) {
  80.           t = wk*gw*(1-gw)/2;
  81.           goto label40;
  82.         }
  83.         else {
  84.           sgn = -sgn;
  85.           wh = wh*wk;
  86.           normbase(&wh, &g2);
  87.           wk = 1/wk;
  88.           if (wk<0)
  89.         b = b+.5;
  90.           b = b-(gw+g2)/2+gw*g2;
  91.         }
  92.       h2 = wh*wh;
  93.       a2 = wk*wk;
  94.       h4 = h2*.5;
  95.       ex = 0;
  96.       if (h4<150.0)
  97.         ex = exp(-h4);
  98.       w2 = h4*ex;
  99.       ap = 1;
  100.       s2 = ap-ex;
  101.       sp = ap;
  102.       s1 = 0;
  103.       sn = s1;
  104.       conex = fabs(con/wk);
  105.       do {
  106.         cn = ap*s2/(sn+sp);
  107.         s1 = s1+cn;
  108.         if (fabs(cn)<=conex)
  109.           break;
  110.         sn = sp;
  111.         sp = sp+1;
  112.         s2 = s2-w2;
  113.         w2 = w2*h4/sp;
  114.         ap = -ap*a2;
  115.       } while (1);
  116.       t = (atan(wk)-wk*s1)/twopi;
  117.     label40:  
  118.       b = b+sgn*t;
  119.     }
  120.     if (is>=0)
  121.       break;
  122.       } while(ak!=0);
  123.     }
  124.     else if (r>=0)
  125.       if (ah>=ak)
  126.     b = 2*gh;
  127.       else
  128.     b = 2*gk;
  129.     else if (ah+ak<0)
  130.       b = 2*(gh+gk)-1;
  131.   }
  132.  label50: 
  133.   if (b<0)
  134.     b = 0;
  135.   if (b>1)
  136.     b = 1;
  137.  
  138.   return(b);
  139. }
  140.